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A new mechanism for relaxing the nuclear magnetic moments, in which a pair of spin-polarized 
BCS quasiparticles, is emitted or absorbed, and which dominates at low temperature, is identified in 
type-II d-wave superconductors in an external magnetic field above H c \ . The results of the theory 
are compared with the NMR experiments on YBCO in high magnetic fields and found to agree 
without invoking antiferromagnetic order in the vortex core. 



I. INTRODUCTION 

Being a bulk real-space probe, with information about 
the precession and relaxation rates of nuclear spins at dif- 
ferent sites inside the sample, nuclear magnetic resonance 
(NMR) has served as one of the key experimental toolsi 
in the study of the electronic properties of high temper- 
ature cuprate superconductors (HTS). By and large, the 
17 NMR data^r— on HTS, in moderate to large magnetic 
fields, has been interpreted as evidence for antiferromag- 
netic (AF) order in the vortex core^— . This interpreta- 
tion was based on the dependence of the spin-lattice re- 
laxation rate, 1/Ti, for the Cu nuclei in the normal state; 
the rate remains constant with temperature^, as opposed 
to being linear in T, as predicted by the Korringa lavA 

We reexamine this interpretation using analytical and 
numerical solutions of the the Bogoliubov-de Genncs 
equations and find that the present data can be under- 
stood quantitatively without invoking AF ordering. In 
particular, the low temperature upturn in 1/(T\T) near 
the vortex core, but not away from it, can be under- 
stood to be caused by a combination of two effects: 1) 
the increase of the quasiparticle (qp) wavefunction near 
the core and 2) the shift of the minimum of the qp band 
with spin along the applied magnetic field (spin up) to 
negative energies due to Zeeman coupling. As a result, a 
new electronic channel opens up for relaxing the nuclear 
spin whereby a pair of spin up quasiparticles is emitted or 
absorbed. This is in contrast to the standard "spin- flip" 
channel, in which a spin up (down) qp is destroyed and 
a spin down(up) qp is created. At temperatures below 
the Zeeman scale the latter channel freezes out since the 
number of qp's with spin anti-aligned with the -B-field be- 
comes thermally activated, and the former channel domi- 
nates. In addition, we find that the broad NMR lineshapc 
appears even in the extreme type-II limit, where the dia- 
magnctic response of the 77-field induced supercurrents 
can be neglected. The broadening here is found to be 
due to the spatially non-uniform paramagnetic response 
of d-wave superconductors in the vortex state. 

Several other theoretical investigations of NMR in the 
mixed state of the cuprate superconductors have been 
carried out. One of the earliest treatments was that of 
Takigawa et. al*£ using a self-consistent method of solu- 



tion of the Bogoliubov-de Gennes equations due to Wang 
and MacDonaldii. They found that 1 jT\ is linear in tem- 
perature near the vortex cores at low temperatures and 
exhibits a small, Hebel-Slichter-like, peak near the super- 
conducting transition temperature. At low temperature, 
the rates near the core are also found to be larger than the 
rates away from it, which approach the usual T 3 depen- 
dence. NMR in the d-wave mixed state was also studied 
using a semiclassical approach^, and using a linearized 
form of the Bogoliubov-de Gennes equations^ 3 -. The re- 
sults of the linearized model give faster rates near the vor- 
tex cores than away from them. They found that 1/(T\T) 
near the core increases slowly with temperature up to 
30 K, and remains almost constant over the same tem- 
perature range away from the core. Importantly, these 
works focused on the quasiparticle spin- flip channel, but, 
as mentioned above and as we show in more detail below, 
the Zeeman coupling of the quasiparticles, which cannot 
be ignored at large magnetic fields, introduces an addi- 
tional channel for spin-lattice relaxation which is found 
to dominate at low temperature. 

Our paper is organized as follows. In Sec. II, we state 
the full Hamiltonian of our system and review the basic 
formulas for the Knight shift and 1/T\T. In Sec. Ill, 
we introduce our model for the electronic contribution to 
the Hamiltonian, derive the formulas for the Knight shift 
and \jT\T for this model, and present the results of our 
calculations. In Sec. IV, we discuss the possible influence 
of antiferromagnetic correlations. Finally, we present our 
conclusions in Sec. V. 



II. THE KNIGHT SHIFT AND SPIN-LATTICE 
RELAXATION RATE 

A. Basic Model 

We will start by stating the Hamiltonian describing 
our system. 

H = H e + H n +H hf , (1) 

where H e is the Hamiltonian of the electrons on their 
own, H n is that of the nuclei on their own, and Hhf 



is the hypcrfinc interaction between the electrons and 
nuclei. The nuclear contribution is just the total energy 
of the nuclear spins in an applied magnetic field, 
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where 7„ is the gyromagnctic ratio of the nuclei and 
I(r) is the spin of the nucleus at r. There are other 
terms present, such as quadrupole terms and interactions 
among the nuclei^^. However, in large magnetic fields, 
which we will be considering here, these terms are small 
compared to the above magnetic term. The quadrupole 
terms lead to uneven splitting of the nuclear energies, and 
the interactions may lead to a slight broadening of the 
resonances of the nuclei^^. We use the model of Shastry, 
Mila, and Ric o 15 ' 16 for the hyperfine interaction, 



hf 



^ 2 E 



A(r - r')I(r) ■ S(r'), 
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where 7 e and y n are the gyromagnetic ratios of an elec- 
tron and a nucleus, respectively, I and S are their respec- 
tive spin angular momenta, and the coefficients A(r) are 
the form factors for the hyperfine interactio n 15 ' 16 . Note 
that we are taking "f e to have a negative value. The elec- 
tronic contribution will be the subject of the next section. 

To compare our results to experiments 5 -^, we will be 
interested in the relaxation rates and the Knight shifts for 
the in-plane 17 atoms in YBCO. For these atoms, we in- 
clude the contributions to the form factor A(r) from both 
nearest- (n.n.) and next-nearest-neighbor (n.n.n.) cop- 
per atoms: 7 e 7„fi. 2 A(r— r') is equal to 2.317x 10~ 7 eV for 
the n.n. Cu atoms, and 5.794 x 10 -8 eV for the n.n.n. Cu 
atom a 17 ' 18 . As discussed by Zha, Barzykin, and Pinesi 8 -, 
this form factor suppresses contributions from the AF 
correlations to the spin-lattice relaxation rate at O-sites 
in the normal state. We expect this suppression to persist 
in the mixed state, as we will argue in Section IV. 



B. The Knight Shift 

The Knight shift is a change (usually an increase) in 
the nuclear resonance frequency induced by the surround- 
ing electrons^. This can be attributed to an effective 
magnetic field produced by the electrons through the hy- 
pcrfinc coupling to the nucleus. Using first-order time- 
independent perturbation theory on this term and taking 
the thermal average of the result, we get 

H hMf = -7J^I(r) • <5B e// (r), (4) 

r 

where 

5B e// (r)^7eft£A(r-r')(s(r')) (5) 



uj = "f n B, is proportional to the applied magnetic field, 
this means that the resonance frequency is shifted by an 
amount j n SB e ff. 

We will use this formula in the next section to deter- 
mine the Knight shift in a c?-wave superconductor in a 
magnetic field. Note that the above formulas, in general, 
allow for a position dependence of the effective magnetic 
field; the Knight shift in our system will, in fact, be po- 
sition dependent. 



C. The Spin-Lattice Relaxation Rate 

As is well known^, the spin-lattice relaxation rate at r 
is given by 



1 _lE mn Wmn(r)(E m -E n f 
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(6) 



where E n is the energy of the nucleus at r in a state n, 
and it is assumed that E n = 0. We will also assume 
that the energies are equally spaced, i.e. E n — E n -\ = 
5E. This is not exactly true, due to, for example, the 
quadrupole term, but, because we are working in high 
magnetic fields, such contributions beyond the magnetic 
energy are small, and we may treat the eigenstates of this 
term alone as almost exact, which is one of the assump- 
tions made in the use of this formula^. 

The function W mn (r) entering Eq.© is the transition 
rate for the z component (along the H-field) of the nuclear 
spin at site r to go from mH to nh. We can find these 
rates using Fermi's Golden Rule, 

W mn (r) = \ ( J2 \{mQ'\V(r) \nQ) * 5{E mQ . - E nQ ) 
\QQ> 

where V(r) is the hyperfine interaction in the form, 

V(r) = - leln h 2 J2 Mr - r')i(r) ■ S(r'). (8) 



We will employ these formulas to determine the spin- 
lattice relaxation rate in a d-wavc superconductor in a 
magnetic field in the next section. Again, note that our 
formulas, in general, allow for a position dependence of 
the relaxation rate and the Knight shift. 



III. NMR IN A d-WAVE SUPERCONDUCTOR 



A. The Bogoliubov-de Gennes (BdG) Equation 



is the effective magnetic field experienced by the nu- 
clei and produced by the electrons; (.) denotes a ther- 
mal average. Because the nuclear resonance frequency, 



We now discuss the electronic contribution to the 
Hamiltonian. Our starting point will be the Hamilto- 
nian for electrons on a square tight-binding lattice in a 



magnetic field with a singlet pairing term 
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where the tight-binding hopping constants, t rr i are t rr / = 

—te^ lA "' , A rr ' = Jl A(r) • dr; A(r) is the vector po- 
tential associated with the (constant) applied magnetic 
field B, \g[iB is the spin magnetic moment of an elec- 
tron, and [i is the chemical potential. For the symmetric 
gauge A(r) = |B x r, the values of A VT < relevant for 
a square lattice are -A r r+ x = — 7ry<j>/<f> and -4 rr+ y = 
7rx$/$o, where $ is the magnetic flux through a plaque- 
tte, and $o = hc/e is the flux quantum. The pairing 
field A rr ' is assumed to have a constant magnitude and 
A rr ' = rjs Aoe* 61 "-"-' , where rj±% = —T]±y = 1, otherwise 



i]s = 0. The phase factor e 
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where 



| e i*(r) +e .*(r')| 

</>(r) satisfies the equations, V x V<fi = 27rz J^i — r i)) 
where the are the positions of the vortex cores, and 
V 2 = 0. These conditions determine </>(r) up to terms 
of the form tfio + \r Q ■ r; these constants are fixed by requir- 
ing zero overall current. The vortex cores form a periodic 
Abrikosov lattice, such that each primitive cell (magnetic 
unit cell) of this lattice carries exactly one quantum of 
magnetic flux hc/e. Note that our assumption of a con- 
stant magnetic field effectively corresponds to an infinite 
penetration depth. We choose to assume a pairing field of 
constant amplitude, placing all of the vortex physics into 
the phase. We do so because we believe that assuming 
a constant amplitude, as opposed to calculating it self- 
consistently, will not greatly affect our results. Lacking a 
microscopic model for cuprate superconductors, it is un- 
certain whether a self-consistent calculation will result in 
much improvement of our results. Finally, for notational 
convenience we denote the Zeeman factor by h = ig/i^-B. 

Our method of solution for this problem follows 
RefsJ^r— . To diagonalize this Hamiltonian, we 
introduce the singular gauge-Bogoliubov de Gcnncs 
transformatior 
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(10) 

which allows us to rewrite the Hamiltonian in terms of 
the quasiparticles in the Bloch basis corresponding to the 
magnetic unit cell £ x x £ y containing a pair of vortices. By 
Bloch's theorem the Nambu spinors, which are eigenfunc- 
tions of the Bogoliubov-de Gennes equation 
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be written as [ukn(r), «kn(r)] 1 



-,ik-r 



[t/k 



can 

T 



(r),Vu„(r)] . 

where C/kn(r) and Vkn(r) are periodic in i x x £ y , n is the 
band index and the crystal momentum k G (— f-, j-] x 

There is one issue introduced by this transformation 
that is worth addressing in detail. As we wind around a 
vortex, (f>(r) increases by 2ir. This means that the phase 



factors in the above gauge transformation only increase 
by 7r, meaning that the phase factors, at this point, are 
not uniquely determined. We must therefore introduce 
branch cuts into the Hamiltonian and choose the val- 
ues of the phase factors carefully. The procedure we use 
in choosing the values of these factors is that used in 
Ref^S. We first choose a branch cut, which can be any 
continuous curve connecting the two vortices inside the 
magnetic unit cell. We then choose one point ro on the 
atomic lattice, and let the phase factor for that site be 
bo = e l ^ r °^ 2 . We now move to a neighboring site r such 
that we do not need to cross the branch cut to reach it. 
Let b be the phase factor for this site. The solution to 
b 2 = e 1 ^*' that we choose is the one that gives the lower 
value of |6 — &o|- We do this for all sites, thus generat- 
ing the appropriate values for e l ^ r ^ 2 . This process is 
illustrated in Figure [U 




FIG. 1: Illustration of the process by which we choose the 
values of e"^ r ' / ' 2 for each site for a 6 x 6 magnetic unit cell. 
The large circles represent the sites, and the small circles the 
bonds connecting them. The blue dots are the locations of 
the vortex cores, and the dashed blue line is the branch cut. 



The black arrows represent the site phase factors, e 



i<t>(r) 



while 



the orange arrows represent the value of e ll?> ^ r '' /2 chosen by 
the process outlined in the text. The bond phase variable, 
e i& "' , is represented by the gray arrow, and e ^( r )/ 2 e ^( r ')/ 2 
is represented by the red arrow. 



The coefficients ?7kn(r) and Vk„(r) satisfy the 
Bogoliubov-de Gennes equation, e _lk l "Hoe lk ' r< f ) kn(i') = 
£kAn(r), where $ k „(r) = [t/ k „(r), Vk„(r)] T , H = 



CTz(£ r — m) + Ca;A r — h, the operators £ r and A r are 

S=±x,±y 

A r = A 2J Z2.r,r+5VsTs, 
<5=±x,±y 

and T,5 performs a translation along the vector 5. 

We now wish to make a comment on the energies of the 
quasiparticlcs in our system. Let Hq^nz = cr z (£ T — fi) + 
a x A T — that is, Hq.nz is Hq without the Zeeman term. 
Note that diagonalizing "Ho.jvz is the same as diagonal- 
izing Jio because the two differ only by a term propor- 
tional to the identity matrix. In fact, if we let Ek„ be the 
eigenvalues of T-Io,nz, then the eigenvalues of Jio are just 
E = Ek n — h. The matrix Hq.nz, as we implied, would 
replace Ho if we neglected the Zeeman splitting. Since 
we can simultaneously diagonalize Ho and Hq,nz, we see 
that the same eigenvectors would diagonalize the differ- 
ence between the two Hamiltonians that result in these 
matrices, which is a term proportional to the z compo- 
nent of the spin. This is because the z component of the 
spin is a good quantum number, and can be used to label 
the elementary excitations. We note that 



[{i<T y yHo,Nz{-iv V )]* =<?y'HQ,NZ a y 



-a z {£ T - fi) - a x A T = -H 



0,NZ- 



(11) 



This shows that, if we multiply an eigenvector of Ho.nz 
by icr y and then take the complex conjugate of the result, 
then we obtain another eigenvector of the same matrix, 
but with the negative of the eigenvalue of the original vec- 
tor. This is exactly what we did to obtain V'kn(r) from 
4j'{y). We have thus shown that the spinor ip'(r) is an 
eigenvector of Ho.nz with eigenvalue — Eknt an d there- 
fore also an eigenvector of Ho with eigenvalue —Ek„ — h. 

The fact that we can generate the negative-energy 
states from the positive-energy ones implies that we may 
take the sum in equation (|10p to be over states that give 
positive eigenvalues of Hq,nz- We could also choose, for 
example, the states with negative eigenvalues. In fact, 
we may choose any one of these two possibilities for each 
term. In this paper, we conform to widely-used conven- 
tion and use the positive energy eigenvalues. 

The diagonalized Hamiltonian takes on the form 

H = X>k„t7 k „ t 7knt + £kni7 kr4 7kni) - £ (0) (12) 

where E^ — Nfi + X^kn an( l the IP eigenenergies 
are E\^ na = E\^ n — ah. The density of states (without 

Zeeman) N(oj) = J2n=l I sTS7^( w — ^kn) for realistic 
values of the physical parameters is plotted in Fig. [5] 
As shown below, these energies and wavefunctions enter 
into the calculation of the NMR line shape and 1/T±. 

In our calculations, we worked with both 20 x 34 and 
26 x 26 unit cells, both of which correspond roughly to 
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FIG. 2: Density of qp states for B = (green line) and B « 
42T (£ x = 20a, £ y = 34a; black line) for a D = 14 and p. = 
0.3t. The vertical line shows the corresponding Zeeman shift 
h = ^gfJsB. The inset shows the spin polarized qp Fermi 
surfaces which come from the lowest (red) and the next to 
lowest (black) magnetic bands. 



an applied field of 42 T, as well as with a 36 x 62 unit 
cell, which corresponds to an applied field of about 13 T. 
These calculations were done for optimally-doped YBCO, 
for which t = 153 mcV, A = ±t, and jU = 0.297£. 
Because we wish to calculate thermodynamic properties 
of the system at low temperatures (temperatures up to 
30 K), we only needed to find some of the lower energy 
bands. Thus, we used the Arnoldi method to find the 
energies and wave functions, and we discretized the re- 
ciprocal lattice into a 50 x 50 grid. 

Let us now make some comments on the energy spec- 
trum and the wave functions. First, we note that, while 
the energies E-^ n are all positive, so that no quasiparti- 
cles would be present in the ground state of our system if 
there is no magnetic field, it is possible, under certain cir- 
cumstances, for some of the energies to be negative. 
To be exact, E^n for most values of the chemical poten- 
tial could exhibit a gap. If £k n t = -Elm — h is smaller 
than the Zeeman splitting, then some of the energies will 
become negative^. This means that, in the ground state, 
there will be some quasiparticlcs present, all with their 
magnetic moments parallel to the field. This means that 
the "gas" of quasiparticlcs is spin-polarized in the ground 
state, resulting in a non-zero Knight shift, even at zero 
temperature. 



B. The Knight Shift and Line Shape Broadening 

To find the Knight shift for a superconductor in a mag- 
netic field, we simply substitute the Bogoliubov transfor- 
mation (jlpp into the spin operator. Here, we only con- 
sider the effect of the spin component along the z axis, 

&(r) = 4(r)St(r)-eJ(r)^(r). 
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FIG. 3: (a) Effective magnetic field shift SB e ff at T = as 
seen by the 17 O nuclear spins for (approximately) triangular 
vortex lattice corresponding to the external field B — 42 T. 
The Dirac anisotropy qd = t/Ao = 14 and fi = 0.3t corre- 
sponding to x ~ 15%. (b) Spatial variation of SB e ff, whose 
density is shown in red, leads to broadening of the NMR line 
shape (black) (additionally broadened by a Gaussian with 
o = 50 gauss^). (c) Spin-lattice relaxation rate 1/(T\T) vs. 
T for different SB e ff. The data points and error bars are the 
experimental data-. 
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FIG. 4: (a) Effective magnetic field shift 5B e ff at T = as 
seen by the 17 nuclear spins for (approximately) triangular 
vortex lattice corresponding to the external field B — 13 T. 
All other parameters are the same as before, (b) Spatial vari- 
ation of SB e ff, whose density is shown in red, leads to broad- 
ening of the NMR line shape (black) (additionally broadened 
by a Gaussian with a — 50 gauss^). (c) Spin-lattice relax- 
ation rate \/{T\T) vs. T for different 8B e ff. Again, the data 
points and error bars are the experimental data-. 



Upon performing the Bogoliubov transformation and 
taking the thermal average, we find that the effective 
magnetic field shift is 

SB ef f (r) = 7e ft££A(r-r> kn (r')[/(£kni) - /(^knf 

r' kn 

(13) 

where f(E) is the usual Fermi-Dirac distribution, 

- cB /Jt + 1 . ( 14 ) 

and n k „(r) = |w kn (r)| 2 + |v k „(r)| 2 . 

Note that, even in the extreme type-II limit (i.e. tak- 
ing the penetration depth to infinity), in the vortex state 
the local electron density is different at different loca- 
tions within a magnetic unit cell. Taking into account 
the Zeeman shift, this translates to spatially varying spin 
density and by the above two Eqs. to the spatially vary- 
ing Knight shift. This means that even nuclei of the 
same species will have different resonance frequencies de- 
pending on their location in the magnetic unit cell. This 
results in a broadening of the NMR line shape (Figures 
12 and @|). 

We found the effective magnetic field due to the elec- 
trons, which is proportional to the Knight shift, for the 
20 x 34 and 36 x 62 lattices, and we have plotted the spa- 
tial profiles for this case in Figures [3] and 01 As expected, 
the largest Knight shifts occur near the vortices because 
the local electron density is largest in the same area. In 
reality, the plots shown are for those O atoms on bonds 
parallel to the "short" axis of the magnetic unit cell (in 



this case, the x axis); the plots for the atoms on bonds 
along the y axis are similar in appearance. 

We now determine the line shape that would result 
from this effective magnetic field. Denoting this line 
shape by f(Bo), the line shape is given by 

f(B ) = J S[B - B(r)} d 3 r. (15) 

This formula assumes that the response of a single nu- 
cleus as a function of the frequency is a delta function 
centered at the resonance frequency. This, however, is 
not true in reality; in fact, the response has a finite width 
in the frequency. For this reason, we must convolute this 
line shape with a broadening function representing the 
response of a single nucleus to obtain the true line shape. 
The broadening function we use in our calculations is a 
Gaussian of width 50 G. This broadening is experimen- 
tally motivated; the width is approximately that of the 
-l/2f> —3/2 transition found by Mitrovioil. Again, we 
calculate this line shape for 20 x 34 and 36 x 62 lattices; 
the results are plotted in Figures [3] and [4] 

C. Nuclear Spin-Lattice Relaxation Rate 

To find the nuclear spin-lattice relaxation rate, 1/Ti, 
we first make a simplification to the hyperfine interaction 
©. Using the identity, 7+S1 + I-S + = 2(I X S X +I y S y ), 
we may write it as 

V(r) = - leln h 2 J2c(r - R)[J+(r)S_(R) + J_(r)S+(R)], 

R 

(16) 



where, for convenience, we have defined C = \A. Note 

that we dropped the I Z S Z term; this term will only con- 
tribute when m = n, and these transition rates, as we 
will see, do not contribute to the relaxation rate. This 
form will be more convenient to work with because the 
S± operators take on simple forms, namely S+ = eje^ 

and S- = cjof. 

We now make an approximation. We assume that the 
nuclear Zeeman energy is much smaller than the elec- 
tronic Zeeman energy, and thus we neglect that contri- 
bution to the total energy of the system. This is a good 
approximation if E m — E n <C fcflT. We then obtain 



where J™"(r) = (m| I+(r) \n), and similarly for 7™ n (r). 



This in turn gives the main result of this paper, 
1 
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The above integrals arc over the 1 st Brillouin zone whose 
area is £Ibz = 47r 2 / '{i x (-y) ■ The qp coherence factors 
enter via the functions G^ k , n , (r) = 



Upon expanding out the expression, 

occurring in equation ([?])■ we obtain four terms; only 
two of these will be non-zero, namely the term involv- 
ing (m| I + \n) (n\ J_ |m) and the term of the same form, and Gknk i n'( r ) 



R 



C r -R [K n {K)Uv n ,(R) + V£ n (R)Vkw(R)] e i(k '- k > R 



but with m and n interchanged, 
equation (J7J), we get 



Substituting this into 
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C r -R [V* n (R)Uv n >(R) - £/ kn (R)Vkw(R)] e 1 
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Note that this differs from the formulas presented in 
Refs i 10 i 13 , most importantly by the presence of the sec- 
ond term. The eigenenergies E nk > are the solu- 
tions of the Bogoliubov-de Gennes equation without the 
+(m <-> n, Q <-> Q))S(Eq/ — Eq ))(18)Zeeman coupling, and the corresponding periodic wave- 



(<m| 7+(r) |n) (n| L(r) |m) (Q'| S_(R) \Q) (Q\ S+(R) |g') 



At this point, we express the spin raising and lowering 
operators in terms of the quasiparticle operators using 
the above definitions and equation (|10[) . and then intro- 
duce these operators into the above expression. Upon 
doing so, we obtain 16 terms. However, 10 of these will 
be zero because they will involve expressions such as 
(Q'\ll\Q) (<3I 77^ \Q)i an d it is impossible to "match" 
the operators in the first factor to those in the sec- 
ond — that is, we cannot pair, for example, an an- 
nihilation operator in the first factor with its corre- 
sponding creation operator in the second. Of the six 
terms that remain, one of them, which has the form, 

(Q'\ lili \Q) (Q\ lilt \Q')i wm a l so be zero because the 
process of creating or destroying two spin down quasipar- 
ticles violates conservation of energy due to all spin down 
quasiparticles having positive energy. We note, however, 
that the corresponding process for spin up quasiparti- 
cles, (Q'\lflf \Q) (Ql^^ \Q'), does not violate conser- 
vation of energy because some of the spin up quasiparti- 
cles have negative energies. This means that, in addition 
to the usual spin- flip scattering process (SF), there is 
also a quasiparticle creation/annihilation process (PC A) 
through which the nuclear spins can relax. 

We go into detail on how we find the different terms 
occurring in our final result in the Appendix; we only 
quote the final result here. Performing the appropriate 
sums and thermal averages we eventually obtain 

W mn (r) = 27r 7e 2 7 ^ 3 [/""(r)/"™(r) + c.c] /(r,T), 



functions are normalized within the magnetic unit cell: 
E re ^ B (l^k„| 2 (r) + |Vk„| 2 (r)) = 1. From Eqs. ®- 
([2"0]) it is readily seen that, regardless of the minimal value 
of £„k, at temperatures T <C h = ^g/isB, the qp spin- 
flip (SF) process is activated and thus vanishingly small. 
At B = 42 T this corresponds to a temperature scale 
of ~ 28 K, which in turn means that the low T(~ 5 K) 
upturn in l/(T\T) observed experimentally 3 ^ cannot be 
due to this process. It is the second term (PCA) which 
dominates at low temperatures and corresponds to the 
observed effect. 

To illustrate the basic physics behind the effect, we 
will temporarily ignore the orbital effects and analyze 
the consequences of the Zeeman coupling alone^. Phys- 
ically, this would correspond to a thin film in a par- 
allel (in-plane) B-field. The eigenenergies in Eq. (f2"0")) 
are then easily found to be E k = \J £^ + A 2 ., where 
6k = — 2i(cos k x a + cos k y a) — /i and Ak = 2Ao(cos k x a — 
cos k y a). At the same time the wavefunctions are simply 



Ik. „ nr i ,„ _ sgnA k 
Ek and «k - ^ 



1 



#- . Assuming 



for simplicity C r = co^ r ,o and particle-hole symmetry, 
this gives for the Zeeman-only case 



1 



— = 27r 7 2 7 2 S 3 c 2 x 



dE 



N(E)N(E + 2h) 



4 cosh 



' E+h ^ 
, 2T i 



21, 



dE 



N(E)N(2h - E) 
8 C osh 2 (^A) 



(21) 



For h,T <C A we need only the low energy qp density 
of states, which is N(E) — 2E/(itvpva), where vp = 

2y/2aU 



'1 — and va = vpAo/t. In this limiting case, 
the integral can be found in a closed form and 



1 



4 llll^cl 



T F 



J F U A 



(22) 



where Fix) 



+ 8xln(l + e x ) - 3a; 2 + &Li 2 (-e x ) 



and Li s (z) is the polylogarithm. For x <C 1, F(x) = 

^ x , and in this limit we recover the standard d- 

wave l/(TiT) ~ T 2 . On the other hand, for x > 1, 
F(x) = a; 2 - z£. In this limit 1/(T X T) increases as T 
is lowered and approaches a constant at T = 0. The 
minimum in 1/{T\T) then results from the competition 
between the spin-flip process which dominates at T > /i 
and the qp pair creation/ annihilation process which dom- 
inates at T <C h. 

Putting back the coupling of the £?-ficld to the orbital 
motion of the electrons, we find that the effect described 
above acquires an interesting spatial content. The dis- 
persing states which are pulled below zero energy by the 
Zeeman coupling are strongly concentrated around the 
cores. Due to the increase of the low-energy wavefunc- 
tions near the cores, the low-temperature relaxation rate 
of the nuclear spin is the largest in the vicinity of the cores 
and decreases with increasing T. This temperature de- 
pendence is in turn due to the pair creation/annihilation 
processes, i.e the second term in Equation (f2"U)) . 

To generate a dependence of the spin-lattice relaxation 
rate on the internal magnetic field, we associated the rate 
at a given point and a given temperature with the value 
of the effective magnetic field shift at that point and tem- 
perature. We then fit this set of points to a power law to 
generate a continuous dependence; at all magnetic field 
shifts of interest, the points are close enough together 
that they approximately form a continuum. We do this 
for all temperatures up to 30 K for B = 42 T and up 
to 20 K for B = 13 T; the results of this procedure are 
shown in Figures [3] and SJ We also highlight the curves 
that give the best fit to the data near the vortex core and 
away from the core^ 



IV. CONTRIBUTION OF 
ANTIFERROMAGNETC CORRELATIONS TO 
THE SPIN-LATTICE RELAXATION RATE 

We now address the issue of how much of an effect an- 
tifcrromagnctic correlations will have on the spin-lattice 
relaxation rate of 17 in the vortex state, assuming that 
the vortex cores represent normal-state regions. As was 
mentioned before, it is known that the form factor filters 
out such correlations in the normal stated. Within the 
simple model presented here, we find that this filtering is 
still active in the vortex state. To investigate the effect of 
vortices, we used a modification of the phcnomenological 
model set forth, among others, by Zha, Barzykin, and 



Pinesi^. We start with their expression for the "antifcr- 
romagnetic" part of the susceptibility, 



XAF(k,w) = 



iy 

4 4^ 1 + (q 



Qi) 2 + iu/usc 



(23) 



Here, a is a scale factor, £ is the antiferromagnetic corre- 
lation length, /ib is the Bohr magneton, the Qi are the 
locations of the peaks in the susceptibility found from 
neutron scattering experiments, wsc is the characteris- 
tic frequency of spin fluctuations, and q ranges over the 
entire first Brillouin zone^. We obtained the model we 
used by rewriting the above susceptibility as a function 
of position, separating the position dependence into a de- 
pendence on the position of a "magnetic unit cell" and a 
dependence on position within the cell, and Fourier trans- 
forming the result with respect to the magnetic unit cell 
position. The result is xaf (q ; Sr — Sr' ,uj) = 



ljV«y ^ _ a£ 2 t i B e l( - cl+G) - {Sr - Sr "> 



+ (q + G - Qi) 2 + jw/wsc ' 



(24) 



where q now ranges over the first magnetic Brillouin 
zone, G is the set of all vectors such that e lG R = 1 
for all R in the magnetic lattice and such that q + G is 
within the first atomic Brillouin zone, Nm is the number 
of magnetic unit cells, and N is the number of sites in 
the atomic lattice. We note that, assuming an L. x x L y 
magnetic unit cell, the number of atomic sites is just 
L x L v Nm, so we may simply write Xaf (Qi Sr — Sr 1 ,oS) = 



1 1 



4 L x L y 



G i 



ae 2 /iae'(q+ G Hfr--fa') 
+ (q + G - Qi) 2 + iuj/u sc ■ 



(25) 



So far, we have not introduced a new model; we only 
rewrote the original in a more complicated form. We 
will now modify this model to introduce vortex effects. 
For simplicity, we will assume only one vortex per unit 
cell. We model these effects by simply introducing step 
functions into the susceptibility to restrict the antiferro- 
magnetic correlations to within a distance equal to the 
superconducting coherence length £sc from the center of 
the vortex. The result is XAF(q, 5r, Sr' , uj) = 

4 L x L y ^ 2f 1 + (q + G 



Qi) 2 + ioj/u>sc 
x0(£sc-<fr)0(£ sc -<5r'). (26) 

We now turn our attention to finding the spin-lattice re- 
laxation rate from this susceptibility. It can be shown 
that, if we assume that this relaxation is due to the hy- 
perfine interaction ((SJ), then the spin-lattice relaxation 
rate will be 

^i- = jhn^ksT J2 A (r- r')A(r -r")x 

Xlim Di *+-(^» ; 



where Di/(. . . ,ui) is a "discontinuity" operator, defined 
as 



Di/(...,w) = 



+ i0+) -/(..., w-iO-) 



2/ 



(28) 



We may rewrite the above formula for a susceptibility 
of the form we are working with. By introducing the 
appropriate Fourier transforms, we eventually arrive at 
the desired result, 



1 



Ti(r) 



7 2 7 2 ft 3 fc B T 



1 



L*L*N M 



xe 



J2 J2 X] A ( k+Gi ) j4 *( k+G2 ) e<(Gl ~ G2) ' r> 

i5r'5r"GiG 2 k 

-i(k+Gi)-<Sr' e i(k+G 3 )-ir" y Di X+- (k, <fr , Sr , u) 



w-»0 



So; 



(29) 



where or' and Sr" are summed over the entire unit cell, k 
is summed over the entire magnetic Brillouin zone, and 
Gi and G2 are summed over the same set of vectors as 
G in Equation (|25jl . 

The rest of our work was done numerically. We 
used the experimentally-determined parameters given 
by Barzykin and Pines for the susceptibility for 



YBa 2 Cu 3 6j 



We first performed a numerical calcu- 



lation of the rates for r = for both copper and oxygen 
without vortices. The form factors we used are 

A Cu (k) = A + 25(cos k x a + cos k y a) (30) 
Ao,x(k) = 2cos\k x a{C\ + 2C2 cos k y a) (31) 
A , y {k) = 2 cos \k v a (d +2C 2 cosk x a), (32) 

where a is the atomic lattice spacing and the parameters, 
A, B, Ci, and C2 are those given by Zha, Barzykin, and 
Pineal. Because of the fact that the momentum-space 
points (k x , k v ) and (k y , k x ) are both present in the sums 
in our formulas, Equations (|3lj) and (f32|) should give the 
same result. We did this using both Equations (|23|) and 
(|25|) to check our formulas. The temperature range we 
examined was 70 K < T < 300 K. We found that the 
contributions to the rates for both Cu and O decreased 
with increasing temperature, and that the rates for O 
were several orders of magnitude smaller than for Cu; 
such suppression of the rate for O compared to Cu has 
been reported before by Mila and Ricei 6 .. We then re- 
peated this calculation, this time including vortices. We 
set the superconducting coherence length £sc = 2a. We 
found that the temperature dependence of both rates was 
qualitatively the same as before, but that the rates were 
actually enhanced; the copper rates by an entire order of 
magnitude and the oxygen rates by a factor of about 4. 
We believe that this is due to the fact that, by impos- 
ing the distance cutoff, we removed contributions to the 
relaxation rate that would have reduced the rate. Based 
on this simplified model, we therefore expect that, not 
only is the filtering effect still present in the mixed state, 
but it is, in fact, enhanced. 

It goes without saying then that this filtering effect 
would therefore make oxygen NMR more sensitive to 



the pair creation and annihilation (PCA) processes than 
copper NMR. This is because the effect of antiferromag- 
netic correlations is much less for oxygen than for copper, 
meaning that the PCA processes will dominate in oxy- 
gen. 



V. CONCLUSIONS 

In this work, we argued that it is possible to explain 
the broadening of the line shape and the upturn in the 
spin-lattice relaxation rate with decreasing temperature 
observed experimentally^ without introducing antifer- 
romagnetic correlations. The line shape broadening can 
be explained, at least in part, by noting that the Knight 
shift varies with position in the lattice in the vortex state. 
This position dependence leads to each nucleus having a 
different resonance frequency, and therefore to a broad- 
ened line shape. The upturn in the relaxation rate can be 
explained as due to a second relaxation process, namely 
creation and annihilation of pairs of spin-up quasipar- 
ticles, that appears when a magnetic field is applied, 
and this process dominates at low temperatures. We do 
not wish to claim that AF correlations do not exist in 
YBCO, only that certain features of the NMR data once 
attributed to such correlations can be explained without 
them; in fact, there is other evidence for the existence 
of such correlations, namely neutron scattering data^. 
As we argued in Section IV, even in the presence of AF 
correlations, the spin-lattice relaxation rates for O will 
not be greatly affected by them due to the form factor. 

Based on the above arguments we expect that, once 
the vortex lattice melts and the system enters a vortex 
liquid phase, the NMR lines sharpen due to motional 
narrowin g 26 ' 27 . At the same time, we expect that the 
spin lattice relaxation rate, 1/Ti, is determined by the 
faster rates and that the low T upturn persists in the 
vortex liquid. 

This picture, and the density of states shown in Fig|2] 
also predict that if an experiment is performed in a clean 
thin film with a well-ordered vortex lattice in which the 
perpendicular component of the _B-ficld is kept fixed, 
while changing the magnitude of the total B, quantum- 
like oscillations in 1/Ti, due to the oscillations of the 
density of states in energy, would be observed. 

While our calculated line shapes have about the same 
width as the experimental shape for the 42 T case, the 
shapes for the 13 T case have different widths. One possi- 
ble contributing factor to this discrepancy is the fact that 
we neglected the variation of the magnetic field and the 
pairing amplitude over a unit cell. We expect the mag- 
netic field to vary more strongly in the 13 T case than in 
the 42 T case because the vortices are further apart in 
the 13 T case. In fact, in the 42 T case, the distance be- 
tween the two vortices in a unit cell is about 10% of the 
penetration depth, while, in the 13 T case, this distance 
is about 23% of the penetration depth. This variation 
will introduce further broadening, which will be greater 



at 13 T than at 42 T, consistent with our findings. 

We also notice that the "tails" on our calculated curves 
are different in length than those of the experimental 
curves. We believe that this, once again, is due to the fact 
that we neglected the variation of the pairing potential 
over a unit cell. In reality, the order parameter should be 
lower in magnitude near the vortex cores because these 
regions are where superconductivity is beginning to break 
down. This means that we expect our calculated line 
shapes to be more accurate in the lower internal field 
regions than in the high internal field regions. 

Finally, we note that the peak in our curve at 42 T 
is split in two, as opposed to the single peak seen in the 
experimental data££. This suggests that there is another 
broadening mechanism at work besides that due to the 
finite width of the normal-state line shape because such 
broadening can wash out the "split" peak so that only a 
single peak appears. One such possibility is the presence 
of impurities. 

We are able to obtain good fits of our calculated tem- 
perature dependence of the spin-lattice relaxation rates 
to the experimental data using the internal magnetic field 
as our only fitting parameter. We note, however, that the 
values of the internal magnetic fields giving us our best- fit 
curves on the line shape do not quite match the experi- 
mental results. In the experiment, the region away from 
the core was in the vicinity of the peak in the line shaped. 
However, the positions of the corresponding theoretical 
curves do not quite fall on the peak; rather, they are away 
from it. It is possible that this discrepancy may be due, 
in part, to our approximations in solving the Bogoliubov- 
de Gennes equation, and that more realistic modeling of 
the vortex core is necessary to account for this. 
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Appendix A: Details of Derivation of Spin-Lattice 
Relaxation Rate 

Here, we go into some more detail about how we found 
the different terms in the spin-lattice relaxation rate. As 
we mentioned earlier, when we rewrite the spin raising 
and lowering operators in terms of the Bogoliubov quasi- 
particles, we obtain 16 terms, though only 6 give non-zero 
contributions. One term we find that contributes to the 
SF process, suppressing factors of u and v that occur, is 

E E E (Q'\tnlW\Q)(Q\trnt%'rn'l\Q') 
QQ' knk'ji' qmq'm' 

xS(E Q , - £ Q QAl) 



We note that the only terms that will give non-zero con- 
tributions are those for which k = q', k' = q, n = mf, 
and n' = m. We may then write 

E E (Q'\tinl%'n't\Q)(Q\lln>t%nl\Q') 
QQ' knk'n' 

xS(E Q ,-E Q ). (A2) 

We now note that the only non-zero matrix elements will 
be those in which the state Q' is obtained from the state 
Q by scattering a particle from the state with crystal 
wave vector k', band index n' , and spin up into the state 
with wave vector k, band index n, and spin down. This 
means that the energy difference between the two states 
is just Eqi — Eq = Ek n i — i?k'?i't = E^n — E^t n i + 2h. 
This energy difference is independent of the exact many- 
particle states Q and Q' , so we may rewrite the sum on 
these states as a trace: 

E E (<3'l7L;7k'r l 't7k'„'t^l < 3') 

Q' knk'n' 

Y.5{E kn -E Un , + 2h) (A3) 

We may now employ the anticommutation relations 
among the quasiparticle operators to rewrite the above 
in terms of number operators: 

E E (<9'l^(i-Wt)|Q') 

Q' knk'n' 

x8{E VLn -Ey L , n ,+2h) (A4) 

Upon taking the thermal average, the trace becomes a 
product of Fermi functions: 

E /(£kni)(l-/(£W t )) 
knk'n' 

x5(E kn -E k , n , +2h) (A5) 

We will also show the calculation for the term resulting 
from the pair creation and annihilation term because it 
will differ slightly from the calculation given above. This 
term (again suppressing actors of u and v) is 

E E E (Q'l7knt7k'n't|Q)(Q|7U7qWtlQ'> 
QQ' knk'n' q?nq'm' 

x8(E Q , - £ Q ()A6) 

In this case, there are two ways to "match" the operators; 
we may either let (k, n) = (q, m) and (k',n') = (q',m') 
or let (k, n) = (q',m') and (k',n') = (q, m). We thus 
obtain 

E E Wl TkntTW't \Q) (Q\ Tk„t7k'n't + 

QQ' knk'n' 

(Q'lWWtlQ) (Qlt'n'ttnflQ')) 

x5(E Q , - E Q ).(A7) 

The only non-zero matrix elements in this case are those 
in which the state Q' is obtained from the state Q by de- 
stroying two quasiparticles with spin up, one with wave 



vector k and band index n and one ith wave vector k' x6(Ek n + i?k'n' — 2ft.). (A8) 

and band index n'. This means that the energy differ- 
ence Eqi — Eq = — -Eknt — ^k'n't = — -^kn — ^k'ra' + 2ft. 

The rest of the derivation proceeds as before, and we 
eventually obtain 

E (i-/(*w))(i-/(£k'„' t )) 

knk'n' 
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